function [powerMean, powerVariance, nSessions ] = ...
    osl_get_source_power_in_mask(variedInput, maskFile, gridStep, sessionInds)
%OSL_GET_SOURCE_POWER_IN_MASK	average source power within a masked region
%
% [POWERMEAN, POWERVARIANCE] = ...
%                 OSL_GET_SOURCE_POWER_IN_MASK(INPUT_FILES, MASK, GRIDSTEP) 
%   returns a vector of the mean power and variance of the power within a
%   masked region specified by nifti file name MASK, for each subject in
%   INPUT_FILES. GRIDSTEP is used to specify the spatial resolution of the
%   input files. The MASK must have the same resolution as GRIDSTEP. 
%
%   INPUT_FILES can be:
%     - a cell array of nifti file names pointing to enveloped data
%       e.g. {'session1_recon_winavHE_delta2s_ss4mm_ds8mm.nii.gz'}
%     - an oat_recon_results structure for a single subject, the output of
%       the first stage of oat. 
%     - a cell array of file names pointing to oat_recon_results structures
%       saved to disk as .mat files (this happens in the first stage of oat). 
%     - an oat structure, after source reconstruction and/or enveloping
%       have been run. The function will process all sessions which are in
%       the results field of the enveloping or source_recon stages. 
%
% [] = OSL_GET_SOURCE_POWER_IN_MASK(INPUT_FILES, MASK, GRIDSTEP, SESSIONINDS) 
%   only uses a subset of INPUT_FILES, as specified by the indices in
%   SESSIONINDS. 
%
% [POWERMEAN, POWERVARIANCE, NSESSIONS] = OSL_GET_SOURCE_POWER_IN_MASK(...)
%   also returns the number of sessions analysed. 


%	Copyright 2014 OHBA
%	This program is free software: you can redistribute it and/or modify
%	it under the terms of the GNU General Public License as published by
%	the Free Software Foundation, either version 3 of the License, or
%	(at your option) any later version.
%	
%	This program is distributed in the hope that it will be useful,
%	but WITHOUT ANY WARRANTY; without even the implied warranty of
%	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%	GNU General Public License for more details.
%	
%	You should have received a copy of the GNU General Public License
%	along with this program.  If not, see <http://www.gnu.org/licenses/>.


%	$LastChangedBy$
%	$Revision$
%	$LastChangedDate$
%	Contact: giles.colclough@eng.ox.ac.uk
%	Originally written on: GLNXA64 by Giles Colclough, 04-Feb-2014 16:32:49

%% Pass in oil matrix, set of oil reconstruction results structures, or enveloped nifti file names
[isEnvelopingDone, nSessions, envelopeNames, ...
    reconResultsNames, isSingleReconResults] = parse_input(variedInput);

% set default session indices to use all files
if nargin < 4 || ~exist('sessionInds', 'var') || isempty(sessionInds),
    sessionInds = 1:nSessions;
    
% catch incorrect input
elseif islogical(sessionInds), 
    error('Expecting linear (not logical) session indices. \n');

elseif nSessions < length(sessionInds), 
    error('You have asked for more sessions than appear to be provided in the input. \n');
    
% set nSessions to be same length as specified sessionInds
else 
    nSessions = length(sessionInds);
end%if

check_spatial_resolutions(maskFile, gridStep);

%% Read in mask
mask = nii_quickread(maskFile, gridStep);

% check the mask looks like a mask
assert(size(mask, 2) == 1 && max(mask) <= 1 && min(mask) >= 0, ...
       [mfilename ':MaskNotRecognised'], ...
       'The mask file doesn''t appear to be a binary mask file. \n');

% in case hasn't been properly thresholded, e.g. after resampling
threshold = 0.5;
mask      = (mask >= threshold); % converts to logical vector

%% Calculate mean and variance of envelope within mask, for each subject

% to get summary statistics for 
for iSession = nSessions:-1:1, % loop backwards to auto-initialise vectors
    
   % envelope if necessary
   if isEnvelopingDone, 
       envelope = nii_quickread(envelopeNames{sessionInds(iSession)}, gridStep);
   
   else % only source recon done
       if isSingleReconResults, % a structure has been passed in
           envelope = calculate_envelope(variedInput);
       else
           tmp = load(reconResultsNames{sessionInds(iSession)}, 'oat_stage_results');
           envelope = calculate_envelope(tmp.oat_stage_results);
           clear tmp
       end%if structure passed in
   end%if enveloping done
   
   % apply mask
   powerTimeCourse = abs(envelope(mask,:));
   clear envelope
   
   % Store mean and variance of Hilbert envelope power over the masked
   % volume. 
   % To do summary statistics properly for each session, we should
   % concatenate the power timecourses one voxel after another, and take
   % the mean and variance of one long timecourse. 
   % In fact, to calculate the mean and variance, it will not matter the
   % order in which we concatenate. 
   powerMean(iSession)     = nanmean(powerTimeCourse(:));
   powerVariance(iSession) = nanstd( powerTimeCourse(:)).^2;
end%loop over sessions
end%osl_get_source_power_in_mask
%%







%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%% SUBFUNCTIONS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%







%% ------------------------------------------------------------------------
function envelope = calculate_envelope(ReconResults)
%CALCULATE_ENVELOPE calculates enveloped data from recon results structure

fprintf('Calculating source reconstructed data using beamformer weights\n');

ABReconResults = AB_get_source_timecourses(ReconResults); 
voxelData      = cat(1, ABReconResults.source_timecourses{:});
clear ABReconResults

fprintf('Finding envelopes of source data. \n');

time = ReconResults.BF.data.D.time;

windowLength = 2.0; %s - Henry's standard value

envelope = envelope_data(voxelData, time, windowLength);
end%calculate_envelope

%% ------------------------------------------------------------------------
function [isEnvelopingDone,  ...
          nSessions,         ...
          envelopeNames,     ...
          reconResultsNames, ...
          isSingleReconResults] = parse_input(in)
%PARSE_INPUT decides whether input is an oil matrix, recon results
%   structures, or list of nifti files. 
%
%   Either the source reconstruction and/or the enveloping must already
%   have been done. 

% initialise all to false
isReconDone       = false;
isReconResults    = false;
% isEnvelopes       = false;
isEnvelopingDone  = false;
nSessions         = 0; 
envelopeNames     = [];
reconResultsNames = [];

% only two valid structure inputs
isOil                = isstruct(in) && isfield(in, 'paradigm') && isfield(in, 'source_recon');
isSingleReconResults = isstruct(in) && isfield(in, 'BF')       && isfield(in, 'osl_version');

% check whether string inputs are valid
isNifti = @(file) exist(file, 'file') && ~isempty(strfind(file, '.nii'));
isMat   = @(file) exist(file, 'file') && strcmpi(file(end-3:end), '.mat');

if isOil, 
    isEnvelopingDone = isfield(in, 'enveloping')                                               && ...
                       isfield(in.enveloping, 'results')                                       && ...
                       isfield(in.enveloping.results, 'source_space_envelopes_results_fnames') && ...
                       ~isempty(in.enveloping.results.source_space_envelopes_results_fnames);

    if isEnvelopingDone,
        envelopeNames = in.enveloping.results.source_space_envelopes_results_fnames;
        envelopeNames = cellfun(@(name) fullfile(in.source_recon.dirname, ...
                                                 in.enveloping.name, ...
                                                 name), ...
                                envelopeNames, ...
                                'UniformOutput', false);
        nSessions     = length(envelopeNames);
        
    else
        isReconDone = isfield(in.source_recon, 'results_fnames') && ~isempty(in.source_recon.results_fnames);
        
        if isReconDone, 
            reconResultsNames = in.source_recon.results_fnames(oil.source_recon.sessions_to_do);
            reconResultsNames = cellfun(@(name) fullfile(in.source_recon.dirname, ...
                                                         name), ...
                                        reconResultsNames, ...
                                        'UniformOutput', false);
            nSessions         = length(reconResultsNames);
            
        else
            error([mfilename ':NoReconPerformed'], ...
                  ['It appears you have input an oil or oat matrix which ', ...
                   'does not point to source reconstructed results. \n', ...
                   'Please perform source reconstruction and try again. \n']);
               
        end%if isReconDone
    end%if isEnvelopingDone
    
elseif iscell(in), 
    if all(cellfun(@(C) ischar(C) && isNifti(C), in)),
%         isEnvelopes      = true;
        isEnvelopingDone = true;
        envelopeNames    = in;
        nSessions        = length(envelopeNames);
        
    elseif all(cellfun(@(C) ischar(C) && isMat(C), in)),
        isReconResults    = true;
        isReconDone       = true;
        reconResultsNames = in;
        nSessions         = length(reconResultsNames);
        
    else
        error([mfilename ':UnrecognisedCellInput'], ...
              ['Your cell array input was not recognised. \n', ...
               'Please provide a cell array of full file paths ', ...
               'to nifti files holding enveloped data or filenames of ', ...
               'saved oat_stage_results. \n']);
        
    end%if input files exist
        
elseif ischar(in), 
    if isNifti(in),
%         isEnvelopes      = true;
        isEnvelopingDone = true;
        envelopeNames    = {in};
        nSessions        = 1;
        
    elseif isMat(in),
        isReconResults    = true;
        isReconDone       = true;
        nSessions         = 1;
        reconResultsNames = {in};
        
    else
        error([mfilename ':FileDoesNotExist'], ...
              'The specified file %s does not exist. Have you used a full path?\n', ...
              in);
    
    end%if input file exists
    
elseif isSingleReconResults,
    nSessions = 1;
    isReconResults = true;
    isReconDone    = true;
    
else
    error([mfilename ':UnrecognisedInput'], ...
          ['Your input was not recognised. \n', ...
           'You may provide an oil or oat matrix after source ', ...
           'reconstruction, a single oat_stage_results structure, \n', ...
           'a list of filenames to oat_stage_results files or a list of ', ...
           'filenames for nifti files holding enveloped data. \n']);
       
end%if is recognised input

if ~nSessions,
    error([mfilename ':NoSessions'], ...
          ['There don''t appear to be any sessions to study. ', ...
          'Something''s gone wrong. \n']);
end%if no sessions

assert(xor(isEnvelopingDone, isempty(envelopeNames)),                             'Inconsistent parsing results. \n');
assert(xor(isReconDone,      isempty(reconResultsNames)) || isSingleReconResults, 'Inconsistent parsing results. \n');
assert(xor(isReconResults,   isempty(reconResultsNames)) || isSingleReconResults, 'Inconsistent parsing results. \n');

assert((isReconDone || isEnvelopingDone), ...
       [mfilename ':InvalidInput'], ...
       ['It appears that your input does not specify enveloped files ', ...
        'or source recon structures. Please check and try again. \n']);
    
end%parse_input

%% ------------------------------------------------------------------------
function isOK = check_spatial_resolutions(mask, res2)
%CHECK_SPATIAL_RESOLUTIONS	Ensures mask and images share the same resolution
[status, res1] = system(sprintf('fslval %s pixdim1', mask));
if status, 
    error([mfilename ':systemCallFailed'], ...
          'System command failed with message: \n   %s \n', res1);
end%if

isOK = isequal(str2double(res1), res2);

if ~isOK, 
    error([mfilename, ':IncompatibleResolutions'], ...
          ['Resolution of mask and images must be the same. \n', ...
           'Try resampling with FLIRT or osl_resample_nii.m. \n']);
end%if
end%check_spatial_resolutions
%% ------------------------------------------------------------------------








% Below: functions by Adam Baker and Giles Colclough, included as
% subfunctions rather than cluttering osl. 

%% ------------------------------------------------------------------------
function [envelopedData, t_ds] = envelope_data(dipoleMag, t, windowLength)
%ENVELOPE_DATA applies Hilbert envelope to data, without normalisation
%
% ENV = ENVELOPE_DATA(DIPOLEMAG, TIME) takes dipoles magnitudes in a set
%    of voxels over TIME, applies a hilbert envelope and downsamples using
%    a window length of 2 seconds. 
%    Dipole magnitudes are expected to occupy rows of DIPOLEMAG: (nVoxels x
%    nSamples). 
%
% ENV = ENVELOPE_DATA(DIPOLEMAG, TIME, WINDOWLENGTH) uses WINDOWLENGTH
%    when downsampling using a moving average. 
%
% [ENV, DS_TIME] = ENVELOPE_DATA(...) returns the downsample time vector
%    DS_TIME of same length as ENV. 


%	$LastChangedBy: GilesColclough $
%	$Revision: 281 $
%	$LastChangedDate: 2014-02-05 12:25:49 +0000 (Wed, 05 Feb 2014) $
%	Contact: giles.colclough 'at' magd.ox.ac.uk
%	Originally written on: GLNXA64 by Giles Colclough, 25-Oct-2013 12:43:45

% Hilbert Envelope
nSamples = size(dipoleMag, 2);
nfft     = 2^nextpow2(nSamples + 1); % zero-padding improves accuracy and speed
HE = abs(hilbert(dipoleMag', nfft)');
HE = HE(:,1:nSamples);

% calculate envelope
if nargin < 3 || ~exist('windowLength', 'var') || isempty(windowLength), 
    windowLength = 2; %s
end%if
overlap          = 0.75;
useHanningWindow = 1;

Fs                  = 1.0 / (t(2) - t(1));
nVoxels             = size(dipoleMag, 1);
winsize             = windowLength * Fs;

for iVoxel = nVoxels:-1:1, % loop backwards to initialise correctly
    % moving average
    envelopedData(iVoxel,:) = osl_movavg(HE(iVoxel,:), ...
                                         t, ...
                                         winsize, ...
                                         overlap, ...
                                         useHanningWindow);
end%loop over parcels
DsFactor = length(t) / size(envelopedData, 2);
newFs    = Fs / DsFactor;
t_ds     = t(1):1.0/newFs:t(end);
t_ds     = t_ds(1:cols(envelopedData));
end%envelope_data
%% ------------------------------------------------------------------------
function results = AB_get_source_timecourses(source_recon_results, varargin)
% Reconstruct source time courses for vector or scalar OSL beamformer
%
% varargin: strings containing additional flags:
%
% 'noreject' - do not reject bad epochs
% 'norecon'  - do not do time course reconstruction
% 'nonorm'   - do not noise normalise
% AB 2013


D=osl_get_oat_sensordata(source_recon_results);

% setup trial list
triallist = D.pickconditions(source_recon_results.source_recon.conditions);

% setuptime indices
t = source_recon_results.BF.data.D.time;
if ~isempty(source_recon_results.source_recon.time_range)
  samples_in_woi = AB_epoch2logical(source_recon_results.source_recon.time_range,t); % fix because woi lumped with bad sections
else
  samples_in_woi = true(size(t));
end
  
if any(ismember(varargin,'noreject'))
  source_recon_time_indices=1:length(source_recon_results.samples2use);
  source_recon_time_indices=source_recon_time_indices(samples_in_woi);
else
  source_recon_time_indices=find(source_recon_results.samples2use);
end


% setup channels
modality = 'MEG';
chanindmeg = strmatch(modality, D.chantype);
chanind = setdiff(chanindmeg, D.badchannels);
if isempty(chanind)
    error(['No good ' modality ' channels were found.']);
end


% load in sensor data
results.sensor_data = D(chanind, source_recon_time_indices,triallist);
results.chanind   = chanind;
results.time_inds = source_recon_time_indices;
results.triallist = triallist;
results.tbad      = AB_epoch2logical(AB_get_bad_sections(D),D.time(samples_in_woi));
results.time      = D.time(samples_in_woi);
results.fsample   = D.fsample;

% HMM stuff
NK=numel(source_recon_results.BF.inverse.W.MEG);
if NK>1
  classchanind=find(strcmp(upper(D.chanlabels),'CLASS'));
  if(isempty(classchanind))
    error(['No ''CLASS'' chanlabel in: ' D.fname]);
  else
    for kk=1:NK
      class_samples_inds{kk} = (D(classchanind,source_recon_time_indices,triallist)==kk);
    end
  end
else
  class_samples_inds{1}= ones(1,length(source_recon_time_indices),length(triallist));
end

results.class_samples_inds = class_samples_inds;


%% Source recon

Nvoxels = numel(source_recon_results.BF.inverse.W.MEG{1});

% Get weights for all voxels & classes:
weights    = cell(Nvoxels,NK);
wnorm      = cell(Nvoxels,NK);
wnorm_nai  = cell(Nvoxels,NK);
variance   = cell(Nvoxels,NK);

for vox = 1:Nvoxels
  for kk = 1:NK
    if strcmp(source_recon_results.recon_method,'none'), % sensor space analysis
      weights{vox,kk}=zeros(1,length(chanind));
      % for sensor space - just compute for all sensors
      weights{vox,kk}(vox)=1;
    else
      weights{vox,kk}=source_recon_results.BF.inverse.W.MEG{kk}{vox};
    end
  end
end


% Get normalisation for all voxels & classes:
for vox = 1:Nvoxels
  for kk = 1:NK
    if ~strcmp(source_recon_results.recon_method,'none') && ~any(ismember(varargin,'nonorm'))
      % this one represents the uncertainty and will be
      % applied equally to the data and regressors, is
      % therefore irrelevant if NK=1
      wnorm{vox,kk} = diag(weights{vox,kk}*source_recon_results.BF.features.C.MEG{kk}*weights{vox,kk}');
      % this one represents a scaling of just the data
      wnorm_nai{vox,kk} = diag(weights{vox,kk}*weights{vox,kk}');      
    else
      wnorm{vox,kk}     = ones(size(weights{1},1),1);
      wnorm_nai{vox,kk} = ones(size(weights{1},1),1);
    end
      variance{vox,kk} = wnorm{vox,kk}./wnorm_nai{vox,kk};
  end
end



results.wnorm = wnorm;
results.wnorm_nai = wnorm_nai;
results.weights = weights;
results.mni_coord = source_recon_results.mni_coord;
results.variance = variance;

if ~any(ismember(varargin,'norecon'))
  results.source_timecourses = AB_get_source_timecourses_recon(results,1:Nvoxels);
end
end%AB_get_source_timecourses
%% ------------------------------------------------------------------------
function dat = AB_get_source_timecourses_recon(recon, voxind)
% AB 2013

% Giles Colclough amended L37: dat{vox,kk}(ori,timeinds,:) -> dat{vox,kk}(ori,timeinds,tri) 

if nargin == 1
  voxind = 1:length(recon.weights);
end

Ntrials = size(recon.sensor_data,3);
Ntpts   = size(recon.sensor_data,2);
Nvoxels = length(voxind);
Nori    = size(recon.weights{1},1);
NK      = numel(recon.class_samples_inds);

dat = {nan(Nori,Ntpts,Ntrials)};
dat = repmat(dat,Nvoxels,NK);

if Nvoxels > 1, ft_progress('init', 'etf'); end

for kk = 1:NK
  
  for tri = 1:Ntrials
    
    % Get subset of data for this this class and trial:
    timeinds = find(recon.class_samples_inds{kk}(1,:, tri)); % time indices for class kk
    if Ntrials == 1 && numel(recon.class_samples_inds) == 1
      sensor_data_subset = recon.sensor_data;
    else
      sensor_data_subset = recon.sensor_data(:,timeinds,tri);
    end
    
    for vox = 1:Nvoxels
      
      if Nvoxels > 1, ft_progress( ((kk-1)*(tri-1)*Nvoxels + vox) / (Nvoxels*Ntrials*NK) ); end
      
      if sum(isnan(squash(recon.weights{kk})))==0,
        if ~isempty(timeinds)
          for ori = 1:Nori
            dat{vox,kk}(ori,timeinds,tri) = recon.weights{voxind(vox),kk}(ori,:)*sensor_data_subset / sqrt(recon.wnorm_nai{voxind(vox),kk}(ori));
          end
        end
      end     
      
    end % vox
    
  end % tri
  
end % kk

if Nvoxels > 1, ft_progress('close'); end

end%AB_get_source_timecourses_recon
%% ------------------------------------------------------------------------
function badsections = AB_get_bad_sections(D,output_format)
% output format: 'logical' or 'indexed'

% AB 2013

% Changed by Giles Colclough 25 Nov 2013 to account for multiple trials. If
% there is more than one trial, D.events is a cell array of structures, one
% for each trial. 

if nargin<2
  output_format = 'indexed';
end

if D.ntrials > 1,
    badsections = [];
    Events = events(D, ':');
    for iTrial = 1:D.ntrials,
        badsections = [badsections; ...
                       get_bad_sections_from_events(Events{iTrial}, ...
                                                    output_format)];       %#ok<AGROW>
    end%if
    
elseif D.ntrials == 1,
    badsections = get_bad_sections_from_events(events(D,1), output_format);
    
else
    error('What''s up with the number of trials?\n');
end%if

end%AB_get_bad_sections



function badsections = get_bad_sections_from_events(Events, output_format)
if ~isempty(Events)
  Events = Events(strcmp({Events.type},'BadEpoch'));
  
  switch output_format
    case 'logical'
      badsections = false(1,D.nsamples);
      for ev = 1:numel(Events)
        badsections = badsections | D.time >= Events(ev).time & D.time < (Events(ev).time+Events(ev).duration);
      end
    case 'indexed'
        duration = [Events.duration]; 
        time     = [Events.time];
        if isempty(duration), duration = zeros(size(time)); end % catch case of duration being empty in each structure (it happens in some of the faces_subject_1 test data!)
        
        badsections = [time; time + duration]';
      
  end%switch

else
  badsections = [];
  
end%if

end%get_bad_sections_from_events
%% ------------------------------------------------------------------------
function l = AB_epoch2logical(ev,t)
% AB 2013

l = false(size(t));

for i=1:size(ev,1)
  l = l | ( t >= ev(i,1) & t < ev(i,2) );
end

end
%% ------------------------------------------------------------------------
% [EOF]